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CONCENTRATIONS OF DARK HALOS FROM THEIR ASSEMBLY HISTORIES 
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ABSTRACT 

We study the relation between the density profiles of dark matter halos and their mass assembly histo- 
ries, using a statistical sample of halos in a high-resolution N-body simulation of the ACDM cosmology. 
For each halo at z = 0, we identify its merger-history tree, and determine concentration parameters 
c v ir for all progenitors, thus providing a structural merger tree for each halo. We fit the mass accre- 
tion histories by a universal function with one parameter, the formation epoch a c , defined when the 
log mass accretion rate dlogM/dloga falls below a critical value 5*. We find that late forming galaxies 
tend to be less concentrated, such that c v i r "observed" at any epoch a is strongly correlated with a c 
via c v i r = c\a /a c . Scatter about this relation is mostly due to measurement errors in c v i r and a c , 
implying that the actual spread in c v i r for halos of a given mass can be mostly attributed to scatter in 
a c . We demonstrate that this relation can also be used to predict the mass and redshift dependence of 
c v i r , and the scatter about the median c v j r (M, z), using accretion histories derived from the Extended 
Press-Schechter (EPS) formalism, after adjusting for a constant offset between the formation times as 
predicted by EPS and as measured in the simulations; this new ingredient can thus be easily incorporated 
into semi-analytic models of galaxy formation. The correlation found between halo concentration and 
mass accretion rate suggests a physical interpretation: for high mass infall rates the central density is 
related to the background density; when the mass infall rate slows, the central density stays approxi- 
mately constant and the halo concentration just grows as i? V ir- Because of the direct connection between 
halo concentration and velocity rotation curves, and because of probable connections between halo mass 
assembly history and star formation history, the tight correlation between these properties provides an 
essential new ingredient for galaxy formation modeling. 



Subject headings: cosmology:theory — galaxies: halos 
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1. INTRODUCTION 

The theory of cold dark matter (CDM) (e.g., Peebles 
1982; Blumenthal et al. 1984; Davis et al. 1985) provides a 
remarkably successful framework for understanding galaxy 
assembly and structure formation in the universe. Within 
this picture, dark matter collapses first into small halos, 
and these halos merge to form progressively larger halos 
over time. As this process continues, baryons that initially 
trace the dark matter cool and condense to form galaxies 
in halo centers. New supplies of gas and galaxy mergers 
are closely linked to the mass accretion histories of the 
halos they inhabit. A detailed understanding of how this 
mass accretion occurs, and how individual halo properties 
depend on their merger histories, is of fundamental im- 
portance for predicting galaxy properties within the CDM 
theory, and, similarly, for using observed galaxy properties 
(e.g., rotation curves) to test the paradigm. 

The basic theory for the buildup of structure in the 
universe, and for the evolution of the properties of 
gravitationally-bound structures, is fairly well developed; 
it has been extensively simulated at increasingly high res- 
olution, and analytic formulations have been developed 
to describe this behavior. The Press-Schechter formalism 



(Press & Schechter 1974; Bond et al. 1991) has provided a 
useful framework for understanding the mass function of 
dark halos, and, with subsequent improvements, has been 
shown to work well in comparison to N-body simulations 
(Gross et al. 1998; Sheth & Tormen 1999; Sheth, Mo, & 
Tormen 2001; Jenkins et al. 2001; Sheth & Tormen 2001). 
The theory has been extended using an excursion-set for- 
malism to an Extended Press-Schechter theory (EPS), 
which describes the statistics of the buildup of individ- 
ual halos through time (Bond et al. 1991; Lacey & Cole 
1993). This EPS theory has been implemented to con- 
struct random realizations of merger trees, each specify- 
ing a full assembly history for a halo (e.g., Lacey & Cole 
1993; Kauffmann et al. 1993; Somerville & Kolatt 1999, 
hereafter SK99). Detailed comparisons of EPS with simu- 
lations have highlighted the general success of the theory 
and quantified the level of accuracy of its predictions (e.g., 
Somerville et al. 2000, hereafter SLKD; Gardner 2001; 
Cohn, Bagla, & White 2001). 

Similar advances have been made in understanding halo 
density structure (Efstathiou et al. 1988; Frenk et al. 1988; 
Dubinski & Carlberg 1991). Navarro, Frenk, & White 
(1995,1996,1997; hereafter NFW), have proposed that halo 
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profiles can be universally fit by a two-parameter func- 
tional form: 

PnfwM = , u f 77TT2' W 
(r/R s ) (1 + r/i? s ) 

where R s is a characteristic "inner" radius, and p s a corre- 
sponding inner density. One of the inner parameters can 
be replaced by a "virial" parameter, either the virial ra- 
dius (i?vir), mass (M v i r ), or velocity (Kir), defined such 
that the mean density inside the virial radius is A v ; r times 
the mean universal density p u at that redshift: 

4-7T 3 

M v i T = — A v i r p u _R v i r . (2) 

The critical overdensity at virialization, A v j r , is motivated 
by the spherical collapse model; it has a value ~ 180 for 
the Einstein-deSitter cosmology, and ~ 340 for the ACDM 
cosmology assumed here. A useful alternative parameter 
for describing the shape of the profile is the concentration 
parameter c v ; r , defined as c v - lr = R v - lr /R s . NFW found that 
this functional form provides a good fit to halos over two 
decades in radius, for a large range of halo masses, and for 
several different cosmological scenarios. They tested it for 
the Einstein-deSitter model with a standard CDM power 
spectrum of initial fluctuations (SCDM), a flat cosmolog- 
ical model with f2 m = 0.3, f^A = 0.7 and a corresponding 
CDM power spectrum (ACDM), and several models with 
power-law power spectra (confirmed by Craig 1997, and 
Kravtsov, Klypin, & Khokhlov 1997). 

Later studies (e.g., Moore et al. 1998; Ghignaet al. 2000; 
Klypin et al. 2001) have indicated that the inner logarith- 
mic slopes of at least some halo density profiles in CDM 
cosmologies are closer to —1.5 than to —1.0. 1 However, 
Klypin et al. (2001) showed that even for halos that are 
better fit by —1.5 slopes, an NFW fit is perfectly adequate 
outside the inner ~ 1% of the virial radius (this corre- 
sponds to roughly <~ 3kpc for a Milky- Way type halo) . An 
important implication (as pointed out by Bullock et al. 
2001, hereafter B01) is that even if halos are not perfectly 
described by the NFW form (see also Avila-Reese et al. 
1999 and Jing 2000), fit parameters derived using this pro- 
file provide a useful general characterization of the density 
structure, relating, for example, the central density of a 
halo to that of the background via two fit parameters (e.g., 
Cyir and _/? v j r ). 

The two-parameter NFW characterization of halo pro- 
files has provided several useful insights into the nature of 
halo collapse. Among the most interesting results, first no- 
ticed by NFW, was that, for a given cosmology, the central 
density p s varies inversely with halo mass, or, equivalently, 
c v j r tends to increase with decreasing M v i r . A natural rea- 
son for this fact is that low-mass halos typically collapse 
earlier, when the universe was denser, and the central den- 
sity somehow reflects this higher background density. 2 In 
a toy model to explain this correlation, NFW (1997) as- 
sumed that p s is a constant multiple k of the universal 
density p u (z c ) at a collapse redshift z c . They defined the 
collapse redshift as the time when half of the halo's mass 
was first in progenitors more massive than / times the 
halo's mass. The general trend of the relation between the 

1 The opposite suggestion has also been made: that the asymptotic sic 
2001). 

2 Indeed, for models with truncated power spectra, in which there is no 
is seen between mass and c v j r (Avila-Rccsc et al. 2001; Bode, Ostriker, 



two profile parameters at z — is reproduced well using 
EPS to predict z c , and a proper choice of values for the 
constants k and / (favored parameters for ACDM were 
/ = 0.01 and k = 3.4 x 10 3 ). 

Subsequent analyses revealed additional complexities 
and trends in halo structure. First, it has been realized 
that while the M v i r -c v i r trend holds on average, halos of 
fixed mass show significant scatter in their c v i r values (Bul- 
lock 1999; Jing 2000; B01). In the context of the proposed 
correlation between concentration and collapse time, the 
scatter in c v i r is a natural reflection of the expected scatter 
in collapse time for halos of a given mass. 

B01 also found that halo concentrations (at fixed mass) 
are systematically lower at higher redshifts, c v j r cx 1/(1 + 
z), implying a much stronger trend than that predicted by 
the NFW toy model. Since by definition the concentration 
at redshift z roughly relates the halo central density to the 
universal background density at that time via c v j r (z) cx 
[p c (z)/ ( o(z)] 1 /3 ; we have c vir (z) cx p c (z) 1 /3/(i + z ). T hc 
observed relation implies that a halo's central density (or 
collapse time) is set only by the halo's mass, independent 
of the redshift when the halo is observed. One reason for 
the shortcoming of the NFW toy model in reproducing 
the proper redshift dependence is that the collapse time 
as defined by this model also depends on the redshift of 
observation. In order to properly match the time evo- 
lution, B01 proposed a modified toy model in which the 
collapse time, denoted now by a c (instead of a c to avoid 
later confusion) is set by the mass only, and is defined as 
the epoch at which the typical collapsing mass, M*(a c ), 
equals a fixed fraction F of the halo mass at epoch a, 
M*(a c ) = FM v i r . The typical collapsing mass is defined 
by a[M*(a)] = 1.686/ D(a), where a(M) is the rms density 
fluctuation on the comoving scale encompassing a mass M, 
extrapolated using linear theory to the present, a = 1, and 
D(a) is the linear growth rate. The implied concentration 
is given by c v i r (a) = Ka/a c , which, with appropriate val- 
ues of F and K, reproduces quite well the dependence of 
concentration on both mass and redshift as measured in 
simulations (sec B01 for details). 

Simplified models of the type discussed so far provide a 
qualitative understanding and useful parameterization of 
the complex processes seen in simulations (for further in- 
sights and an exploration of non- hierarchical power spectra 
see Eke et al. 2001). However, our understanding remains 
somewhat schematic and several important questions re- 
main open. First, how do individual halo density profiles 
build up over time? How do the mass accretion histories 
affect the final halo concentrations, and how can physi- 
cally realistic mass accretion histories be connected to the 
simplified definition of formation time advocated by B01? 
Can EPS be used to predict halo concentrations? What 
physical process is responsible for the scatter in c v ; r at 
fixed mass? This work builds on the work of B01, us- 
ing thc same simulations analyzed there, together with a 
new "structural merger tree" described below. The goal 
is to see if we can directly correlate the assembly history 
of halos with their dark matter halo density profiles, test 

pe as r — > may be somewhat shallower than —1.0 (Taylor & Navarro 

systematic relation between collapse time and halo mass, no correlation 
& Turok 2001; Eke, Navarro, & Steinmetz 2001). 
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the model proposed by B01, and develop a physical model 
for the range of halo concentrations seen in simulations, 
including scatter and dependence on mass and redshift. 

These questions are interesting from a theoretical per- 
spective, but they also have direct observational implica- 
tions. The shape of the halo density profile directly affects 
the rotation curve of the galaxy that forms within it; the 
one-sigma scatter in concentration observed by B01 cor- 
responds to, for a 1 x 1O 11 /i _1 M halo, a scatter in V max 
values of ~ 85 — 105 km/s. Thus scatter in halo density 
profiles is directly related to scatter in the Tully-Fisher re- 
lation, and also has implications for the relative contribu- 
tions of halo and disk to velocity rotation curves. In addi- 
tion, density profiles may affect other aspects of the galaxy 
formation process, such as the efficiency of gas cooling or 
the star formation rate. If halo concentrations are related 
to their mass assembly histories, it may indicate that halo 
concentration is correlated with galaxy type. This would 
have implications for both the zero-point and scatter in the 
Tully-Fisher relation, and could also have implications for 
the rotation curves of low-surface brightness galaxies. Pos- 
sible correlations between halo density profiles and galaxy 
observables (with the exception of mass) have been ne- 
glected in previous modeling efforts, but it is clear that 
such correlations are likely to be quite important. 

We begin in §2 by summarizing the relevant details of 
the N-body simulation, halo finder and density profile fit- 
ting. In §3, the "structural merger tree" developed for 
this project is described. We continue in §4 by detail- 
ing how mass accretion histories are constructed and then 
describe a new method for defining a characteristic for- 
mation epoch for each halo. We then show how this for- 
mation epoch can be related to the halo concentration, 
and how this can explain the dependence of concentration 
on mass and redshift as well as explaining the origin and 
magnitude of the scatter in these relations. In §8, we show 
how Extended Press-Schechter theory can be used to pre- 
dict concentrations for individual halos using our model 
for relating halo concentration to characteristic formation 
epoch. This model can reproduce the scatter, mass and 
redshift trends observed in N-body simulations. In §9, we 
discuss the scatter in c vlT (M v i T ), and how it depends on 
the merging history of halos. We discuss the implications 
of our results and conclude in §10. 

2. SIMULATED HALOS 

In the work that follows, we consider only one cosmol- 
ogy whose evolution has been simulated with the ART 
code (Kravtsov, Klypin, & Khokhlov 1997), a flat ACDM 
model with £l m = 0.3, h = 0.7 and as — 1.0. The sim- 
ulation is the same as that used in B01. It followed the 
trajectories of 256 3 cold dark matter particles within a 
cubic, periodic box of comoving size 60/i _1 Mpc from red- 
shift z — 40 to the present. A 512 3 uniform grid is used, 
with up to six refinement levels in the regions of highest 
density, implying a dynamic range of 32, 768. The for- 
mal resolution of the simulation is thus / res = 1.8h~ kpc , 
and the mass per particle is m p = 1.1 x W 9 h~ 1 M . For 
the purpose of constructing accurate merger trees, here we 
analyze the simulation data from 36 output times thinly 

3 For flat cosmologies, A v i r can be approximated by (Bryan & Norman 
is the ratio of mean matter density to critical density at redshift z; for 



spaced between z = 7 and 0. It should be noted that the 
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Fig. 1. — Output times of the simulation. The curve is the age 
of the universe t as a function of the universal expansion factor a 
for the ACDM model considered here. The symbols mark the 36 
output times of the simulation. 

methods described here for finding and fitting halos and 
constructing merger trees are completely generalizable to 
other simulations or cosmologies. 

The halo finding algorithm used here is based on the 
Bound Density Maxima (Klypin & Holtzman 1997) tech- 
nique described in Bullock (1999) and B01, but we have 
modified and optimized it for the purpose of building a 
structural merger tree. The essential elements are pre- 
sented briefly here, and details are described in the Ap- 
pendix. For each density maximum, we step out in radial 
shells until the mean overdensity falls below A v i r (z), 3 or 
the radial profile shows a significant upturn. We denote 
this radius as i?h and the mass determined by counting- 
particles within this radius as M^. We attempt to iden- 
tify halos containing as few as A™ 111 = 20 particles; our 
halo catalog thus includes ~ 14000 halos above a mass 
threshold of 2.2 x 10 /i M Q . By comparing our sim- 
ulation with a smaller, higher resolution realization of 
the same cosmology, at this mass we estimate our com- 
pleteness to be 70%, and we are ~ 100% complete at 
6.6 x 10 10 h~ 1 M Q . NFW profiles are fit to halos with more 
than Ap™" = 200 particles, corresponding to halos more 
massive than 2.2 x 10 n /i _1 M G . A profile fit to a halo of 
only a few hundred particles carries large errors, but as 
long as the fit is done properly to include these errors, 
reliable concentrations estimates can be obtained (B01). 
Halos in the mass range 2 x 10 11 — 5 x 1O 11 /i _1 M have 
typical fit mass errors of about 10%, and c V i r errors of 15- 
20%. However, a few percent of the time the errors are 
significantly larger than that. We therefore make a rig- 
orous attempt to estimate the errors and take them into 
account in every step of the process. Poor fits are marked 

1998) A vir ~ (18tt 2 + 82a: - 39x 2 )/n(z), where x = Q(z) - 1, and Q(z) 
the ACDM model considered here, A v i r ~ 337. 
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by large errors that are incorporated in the analysis and 
the results we present. The outcome at any output time is 
a statistical halo catalog that includes all the bound viri- 
alized systems in the simulation above the mass threshold 
of 2.2 x lO lo h -1 M . At z = there are 14,219 such halos. 
The output for each halo includes the list of its bound par- 
ticles, the location and velocity of its center, and the NFW 
profile parameters (e.g., c v ; r and M v ; r ) and corresponding 
errors (cr c and um) when applicable. (We also include in- 
formation about the halo angular momentum properties; 
this is not used in the present work, but see also Wechsler 
2001, Vitvitska et al. 2001 and Wechsler et al. in prep). 
The mass function of this revised halo catalog, and a de- 
tailed comparison with the halo catalog used in the work 
of B01, has been presented in Wechsler (2001). 

3. CONSTRUCTING A MERGER TREE 

As a base for the structural merger tree, we use the dis- 
tinct halo catalogs described above at each of 36 output 
times of the simulations, from z = 7 to the present. The 
output times are well spaced in redshift; the cosmological 
scale factors associated with the saved output times are 
shown in Figure 1 . Between each set of output times, we 
track the trajectories of all of the particles. Given a halo 
and an output time, we tag all of its particles and track 
them back to the previous output time. We then make 
a list of all halos at that earlier output time containing 
tagged particles, recording the number of tagged particles 
contained in each one. A similar list of halos and tagged 
particle numbers is obtained for the neighboring future 
output time. In addition, we record the number of parti- 
cles that are not in any halo in the previous output time, 
and the number of particles that do not end up in a halo 
in the subsequent output time. 

We have two criteria for halo 1 at one output time to 
be labeled a "progenitor" of halo 2 at the subsequent out- 
put time. In our language, halo 2 will then be labeled an 
"offspring" of halo 1. First, more than half the particles 
in halo 1 must end up in halo 2. In addition, since, on 
occasion, a halo's particles do not end up in any halo in 
the subsequent output time, we also demand that more 
than 70% of the particles in halo 1 which end up in any 
halo must end up in halo 2. Thus a halo is allowed to 
have only one offspring, but there is no limit on the num- 
ber of progenitors a halo may have. In much of the work 
that follows, we will focus on a "most massive progeni- 
tor" halo in the previous output time. In general, this is 
the progenitor halo which contributes the largest number 
of particles; however if the halo's progenitor mass is not 
at least half of its mass, we additionally require that the 
progenitor's most bound particle is included in the halo. 
Even if this is the case, we also check that the halo's most 
bound particle in the present output time was a member 
of the progenitor; if it was not, the halo is only designated 
as the most massive progenitor if it is more massive than 
the minimum mass of the halo catalog and it is least a 
third the mass of the offspring halo. These criteria were 
designed to maximize the redshift extent of as many mass 
accretion trajectories as possible, without counting trajec- 
tories which may have been affected by the completeness 
of the catalog or severe errors in fitting. 

We have used the criteria outlined above to construct 



the merger tree of every halo at every output time. Ex- 
amples of such a merger tree are shown in Figure 2, which 
shows the mass accumulation history of a small cluster 
halo of mass M vir = 2.8 x 10 14 /i _1 M o at z = 0, and a 
galaxy-mass halo (M vir = 2.9 x 10 12 /i _1 M Q at z = 0). 
Each halo in the tree is represented by two circles, one 
with radius proportional to the halo's virial radius and 
one with radius proportional to the halo's best-fit NFW 
R s parameter. 

4. MASS GROWTH CURVES 

For the purposes of understanding the evolution of halo 
concentrations, it is useful to characterize the history of 
mass assembly in each halo by tracking the evolution of 
its most massive progenitor. The mass assigned to the 
most massive progenitor at each output time is typically 
the best fit virial mass, M v i r . However, for cases in which 
the fit errors on M v ; r were large, we used an iterative proce- 
dure, described in the appendix, to determine a corrected 
mass; this mass is based on either the fit mass, the mea- 
sured mass within i? v i r , or an interpolated mass between 
the previous and subsequent timestcps. In addition, our 
halo mass trajectories do not always extend as far back as 
the first analyzed output time of the simulation. This usu- 
ally happens because the most massive progenitor at some 
redshift falls below our completeness limit and cannot be 
identified, although there are also rare cases in which our 
criteria for progenitor are simply not met (see Wechsler 
2001). In order to have complete trajectories for a reason- 
able number of halos, and in order to have reliable fits for 
most halos, we limit our analysis to halos more massive 
than 10 12 h~ 1 M Q at z = 0. In this mass range, which in- 
cludes <~ 900 halos, <~ 95% of the halo trajectories extend 
back to z = 2, and <~ 90% extend back to z = 3. 

Figure 3 shows the history of mass growth for the major 
progenitors of several different halos, spanning a range of 
masses and concentration parameters. Notice that more 
massive halos tend show substantial mass accumulation up 
to late times, but the growth curves for less massive halos 
tend to flatten out earlier. This tendency is illustrated in 
Figure 4, which shows the average mass accretion histo- 
ries for halos binned by final mass. Again, the high mass 
halos form later than low mass halos, as expected in CDM 
models. 

Since mass accretion is a continuous process, the loose 
term "formation time" is ambiguous and it requires an 
agreed, measurable, quantitative definition. The trends 
of formation time with mass and redshift may depend on 
this definition. One common choice is to set the formation 
time equal to the time when the mass in progenitor halos 
(or in the most massive progenitor) is equal to some frac- 
tion of the halo's final mass M Q (e.g., NFW 1997, Lacey 
& Cole 1993). Definitions of this type have a common 
feature: the formation time is a relative measure which 
depends, for a given halo trajectory, on the redshift z a at 
which the halo is observed. As z a is increased, the forma- 
tion redshift Zf also increases. However, it increases more 
slowly, since mass accretion proceeds more rapidly at high 
redshift in CDM models. Thus both the formation red- 
shift Zf and the ratio of (1 + Zf)/(1 + z ) change with time 
in such a model. As mentioned in the introduction, this 
feature of the formation time definition is the reason that 
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Fig. 2. — Structural merger trees for two halos. This diagram illustrates the merging history of a cluster-mass halo (left; M v i r = 
2.8 X lO 14 h _1 M and c vir = 5.9) and a galaxy-mass halo (right; M vir = 2.9 X 10 12 /i _1 M Q and c vir = 12.5) at a = 1. The radii of the outer 
and inner (filled) circles are proportional to the virial and inner NFW radii, R v j r and R B , respectively, scaled such that the two halos have 
equal sizes at a = 1. Lines connect halos with their progenitor halos. All progenitors with profile fits (M > 2.2 X lO 11 /! - 1 M Q ) are shown for 
the cluster-mass halo; all progenitors (M > 2.2 X 10 10 /i — 1 M„) are shown for the galaxy-mass halo The scalefactor a at the output time is 
listed in the center of the plot. The width of the diagram is arbitrary. 
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Fig. 3. — Selected mass accretion trajectories, showing the evolution of the most massive progenitor for individual halos in the simulation 
(thick). Functional fits to the growth curve of each halo using Eq. 3 are shown as thin smooth lines. The dotted lines represent the expected 
halo trajectory based on the value of the concentration parameter, using Eq. 3, with the value of a c as derived by Eq. 7 using the value c v ir 
measured at z = (and quoted in each panel). The quoted c tra j is derived by Eq. 7 using the a c of the best fit to the actual growth curve. 
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Fig. 4. — Average mass accretion histories, normalized at a = 1. Left: binned in 3 bins by final halo mass: Mq = 1 — 4 X 10 12 h -1 M Q 
(dot-dashed), Mo = .4 — 3 X 10 13 /i -1 M Q (solid), and Mq > 3 x W 13 h~ 1 M Q (dashed). The three (green) curves connect the averages of 
M(o)/Mq at each output time. The pair of dotted lines shows the 68% spread about the middle case (the spread is comparable for the other 
bins). We see that massive halos tend to form later than lower mass halos, whose mass accretion rate peaks at an earlier time. Right: 
binned in 3 bins by formation epoch a c . Dot-dashed lines correspond to early formers (typically low mass halos), dashed lines to late formers 
(typically higher mass halos). The averages and spread are displayed in analogy to the left panel. 



the redshift dependence proposed for the NFW model fails 
to match accurately the dependence seen in simulations. 
Another shortcoming of this kind of definition is that it 
is limited to using the value of the halo trajectory at one 
time, which may introduce noise and miss relevant infor- 
mation. It would be useful to find a quantity which more 
fully characterizes the whole assembly history of the halo, 
and preferably one which does not depend on the redshift 
of observation z Q (at a possible cost that such a definition 
may be allowed to have formation times in the future). 

By examining a range of full mass assembly histories for 
our sample of halos, we have found a useful parameterized 
form that captures many essential aspects of halo growth 
over time. Remarkably, we find that both average mass 
accretion histories and mass accretion histories for indi- 
vidual halos, as observed at z — 0, can be characterized 
by a simple function: 

M(a) = M Q e~ az , a = (1 + z) -1 . (3) 

Although individual halo trajectories may deviate from 
this form significantly in places (e.g., at the time of a ma- 
jor merger), this one-parameter model (in addition to the 
halos' final mass M Q ) provides a remarkably good charac- 
terization of the range of halo mass accretion trajectories. 
Fits to this equation are shown in Figure 3 for several 
representative individual halos. van den Bosch (2001) has 
independently shown that a similar, two-parameter, func- 
tional form can be used to represent halo mass accretion 
histories for a variety of cosmologies and over a large mass 
range. 

The single free parameter in the model, a, can be related 
to a characteristic epoch for formation, a c , defined as the 

4 However, S=2 or something similar is required to match the B01 model 



expansion scale factor a when the logarithmic slope of the 
accretion rate, dlogM/dloga, falls below some specified 
value, S. The functional form defined in Eq. 3 implies 
a c = a/S. 

The same formation epoch can be defined equivalently 
for any "observing" epoch z Q of that halo, by replacing a 
in Eq. 3 by a/a a . In this case, the characteristic formation 
time is related to a via 



a a a I S. 



(4) 



Thus at any such observing redshift, with scalefactor 
a = 1/(1 + z a ) and mass M Q = M(z ), the mass growth 
is fit by 



M(a) = M Q exp 



-a c b I — 
a 



1 



(5) 



This implies that, for any halo whose mass accretion tra- 
jectory resembles Eq. 3, the characteristic formation time 
is the same regardless of the redshift z Q at which the halo 
is observed. In what follows we have chosen S = 2. Since 
the value of S is not used for the fit but only serves to 
define a c , this choice is arbitrary. 

The range of formation scalefactors defined in this man- 
ner exhibits a log-normal distribution, whose mean value 
increases with increasing mass — low mass halos typically 
accrete their mass early, while high mass halos are typ- 
ically still accreting mass rapidly in the present epoch. 
However, for a given mass, there is a large scatter in for- 
mation epoch. The scatter in the mass accretion trajecto- 
ries for a given mass can be seen in Figure 4; we resume a 
detailed discussion of the mass dependence of a c and the 
scatter about this relation in §7. 

behavior that we describe in 57. 



8 



WECHSLER et al. 




5. CONCENTRATION AND ASSEMBLY HISTORY 

We find that the concentration of a halo is tightly cor- 
related with the characteristic formation epoch as defined 
in the above section. Figure 5 shows the average evolu- 
tion of the concentration parameter for halos in different 
mass ranges, corresponding to the average mass trajecto- 
ries shown in Figure 4. From this figure it is clear that 
halo concentrations have a stronger trend with less scat- 
ter when binned on a c (right) than when binned on mass 
(left). 5 We therefore investigate how c V i r is related to a c 
directly. 

Figure 6 shows the relation between concentration and 
a c for halos at z = 0. The concentration of a halo is 
strongly correlated with its characteristic formation time, 
and a good fit is obtained with the inverse relation: 

Cvir = Ci/dc, (6) 

where c\ =4.1 is the typical concentration of halos form- 
ing today. The scatter about this relation is already not 
too large for all the halos, but we note that most outliers 
fall in one of the following three special categories: 

1. the halo has a truncated trajectory that does not 
extend far back to the past, and thus a c is not well 
determined; 

2. the halo has a significant discontinuity in its 
c V j r trajectory at the final output time only, so 
that this value is not representative of the whole 
trajectory (this can occur if there is a merger 
or other disruption occurring at the final output 
time); or 

3. the assembly history includes a merger that is 
substantially larger than the average accretion rate 



5 Note that the figure only shows this directly for 2 = 0, although it is true at any redshift z when a 
about the average trajectory increases with z, since a halo can't uniquely predict its future. 



at that time, and thus Eq. 3 does not provide a 
good description of the actual history. 

To deal with special case (1), halos with trajectories which 
do not extend as far back as z — 1 are excluded from 
further analysis (fewer than 5% of cases; these halos are 
indicated by plus symbols in Figure 6). For case (2), out- 
lined by squares in Figure 6, we find that a much better 
agreement with the median relation is obtained when the 
last discrepant value of c v i r is replaced by the value of c v i r 
in the preceding output time. We do not attempt to cure 
the problem associated with case (3), except for keeping 
in mind that the outliers remaining in the c v i r -a c relation 
are often due to a failure of Eq. 3 to adequately model the 
history of that halo. With these modifications, the scatter 
in Cvir for a given a c is A(logc V i r ) ~ 0.09, without remov- 
ing additional scatter due to large NFW fit errors (see also 
Figure 9), and A(logc V i r ) ~ 0.05 — 06 when errors in c V i r 
are corrected for. 

We have also examined the dependence of c V h- on the 
merger history of halos, which is correlated with but dis- 
tinct from the mass accretion history of the most massive 
progenitor discussed above. Since halos that did not un- 
dergo a recent major merger are more likely to have ac- 
creted most of their mass early, they typically have earlier 
formation times and higher concentration values (see Fig- 
ure 17). However, we find that the parameter a c we have 
defined based on the mass assembly history is more useful; 
for halos with a given a c value, the occurrence of a recent 
merger is not an important factor affecting the concentra- 
tion. This can be seen in Figure 6, which demonstrates 
that halos which have experienced recent major mergers 
(circles; defined here as a merger ratio of greater than 1/3) 

is measured at z . However, the scatter 



CONCENTRATION AND HALO ASSEMBLY HISTORY 



9 




0.1 0.4 1.0 



Fig. 6. — Concentration versus formation epoch for halos more 
massive than 10 12 h~ 'Mg at z = 0. Halos that had major mergers 
(M2 I Mi/3) since z = 1 are indicated with filled circles, halos that 
did not have major mergers since z = 1 are shown with open tri- 
angles. The size of the points is inversely proportional to the error 
on the points (in both c vlI and a c ). Halos with short trajectories 
are marked with pluses and excluded from further analysis; halos 
whose final c v ir value jumps are outlines with squares and treated 
as described in the text. 
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Fig. 7. — Concentration versus scaled formation epoch a c /a , for 
halos at z = (triangles), 2 = 1 (circles), and z = 2 (squares). The 
4 panels correspond to different mass ranges. At all masses and red- 
shifts, the concentration parameter c v ; r is well fit by the functional 
form c v ; r = cia c /a , where c\ ~ 4.1 (represented by the solid line 
in each panel). 



form later on average but follow the same trend in c V i r vs. 
a c (albeit with more scatter; see §9) as halos with no recent 
major mergers (triangles). 

6. DEPENDENCE ON REDSHIFT OF MEASUREMENT 

We now investigate how the correlation found in the 
previous section between c v i r and a c as measured at z = 
behaves as a function of the redshift z Q . We found that 
c v j r is inversely proportional to a c , or equivalently to the 
fit parameter a from Eq. 3. If the determining factor in a 
halo's concentration is the shape of its mass growth curve 
then we expect that a similar dependence would hold for 
all redshifts. Eq. 4 would then imply that the concen- 
tration should be inversely proportional to a c /a for any 
a a when the halo is observed. In Figure 7, the measured 
concentration values are plotted vs. this scaled formation 
epoch for halos at z = 0, 1, and 2. Halos appear to follow 
the same trend regardless of mass or redshift. 

In order to obtain the most reliable estimate of the pro- 
portionality constant c\ in the linear regression of c v i r and 
a c /a , we use halos in the mass range 1 — 5 x 10 12 h~ 1 M & , 
and consider the errors in both a c and c V i r while exclud- 
ing outlying points that are more than 2-a away from the 
best fit line. Considering halos at z = 0, 1, and 2 together 
yields the same result as obtained for z a halos alone, prop- 
erly scaled by a Q : 

CVir = cia Q /a Cl (7) 

where c\ = 8.2/ S (here we have used S = 2 — ► ci = 4.1.) 
As before, the parameter c\ is the typical concentration 



of halos whose formation time is at the time of measure- 
ment, a c = a D . For the ACDM cosmology considered here, 
and a a = 1, this is the typical concentration of halos of 
~ 7 x 10 13 /i _1 M Q . Figure 7 shows that this formula pro- 
vides a good description of the observed correlation be- 
tween concentration and formation time for halos at all 
masses and redshifts. 

7. PROPERTIES OF THE FORMATION TIME 

As discussed above, B01 showed that the average con- 
centration value measured for halos of a fixed mass scales 
as Cvir oc a . This trend was understood using a sim- 
ple model in which the central densities of halos are set 
by the density of the universe at a characteristic collapse 
time, a c , implying c v - n - oc a Q /a c . Crucial to the success 
of this model was that the definition of collapse time, a c , 
was independent of a for fixed mass halos. In the previ- 
ous sections we showed that a similar result obtains using 
a different definition of the characteristic formation time: 
c v ; r oc a Q /a c , where a c is derived based on the actual his- 
tory of mass accretion in individual halos, rather than on 
a simple universal scaling argument as in B01. In order 
to obtain a consistent trend with redshift, we expect that 
a c (like a c ) will be independent of the epoch z when the 
measurement is performed for a given halo mass. This 
assumption is tested in Figure 8, which shows the depen- 
dence of a c on mass, for halos identified at three distinct 
redshifts in the simulation. Within the errors, it shows 
roughly the same mass trend regardless of redshift; this is 
a key feature that our a c parameter has in common with 
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Fig. 8. — Mass dependence of a c , shown for halos at redshifts z = 0, 1, and 2. Dots represent individual halos. The middle thin lines with 
Poisson error bars indicate median values for the mass bins, and the outer thin lines indicate the 68 % scatter about the median. In each 
panel, the thick solid line is a c as a function of mass from the model of B01 and the thick dashed lines indicate the scatter about this value 
needed to account for the 35% measured scatter in concentration. 



the collapse epoch a c defined by B01. In fact, these pa- 
rameters can be directly associated (for an appropriately 
chosen value of S). For example, we find that a c follows 
the same mass trend as our derived a c (with S = 2) if it 
is defined as in B01 (M*(a c ) = FM) with F = 0.015. 

The value of a c for a given mass halo also shows signif- 
icant scatter, as can be seen in Figure 8. The amount of 
scatter (~ 0.13 in the log) can almost completely account 
for the scatter seen in the c v ; r vs. M relation (discussed 
in detail in §9). Figure 9 shows probability distributions 
of a c and c V i r for a given mass range. The intrinsic scatter 
in c v ir for a given a c is relatively small compared to the 
scatter in c V i r for a given mass, and we find that the mea- 
sured distribution of c v ; r can be practically reproduced if 
Eq. 7 is used to transform each measured a c value into a 
concentration. Thus the scatter in concentration can be 
understood as deriving almost exclusively from the range 
of accretion histories for a given halo mass. 

It is encouraging that our definition of formation time is 
robust both in terms of measuring the same value at differ- 
ent epochs along the growth curve of a given halo, and in 
terms of its average value, for halos of a given mass when 
measured at different redshifts. We explore this further 
in the next section. Although a given halo mass accretion 
history that is a perfect fit to Eq. 3 will have a constant 
value of a c regardless of when it is observed, we find that 
the early part of a halo's trajectory is not a good indica- 
tor of the latter part. This is shown in Figure 10, which 
demonstrates that the formation time measured using the 
first half of a halo's history is not correlated with the for- 
mation time measured using the second half. 

8. REPRODUCING THE RESULTS WITH EPS 



We have demonstrated that the tight correlation found 
in the simulations between a c and c V h- can account for 
the mass and redshift trends and the scatter in c^- u (M). 
However, it would be useful to have a way to model the 
concentrations without such a computationally expensive 
simulation, for incorporation into analytic or semi-analytic 
models. Of course, the model of B01 did this to some ex- 
tent for average halo properties as a function of mass and 
redshift, but this neglects the tight correlation found here 
between individual halo density profiles and their mass ac- 
cretion histories. In this section we compare the mass ac- 
cretion trajectories from our simulation with trajectories 
generated using the EPS formalism and test whether the 
correlation we found, c vn - oc Ac" 1 , can be used to predict 
concentrations for individual halos semi-analytically. 

Comparisons of specific aspects of mass accretion his- 
tories as derived semi-analytically and as extracted from 
simulations have been performed (SLKD; Gardner 2001), 
but here we focus on the quantity relevant for our model- 
ing of Cvir, namely the range of values of a c (M, z), which 
has not been investigated previously. There are essentially 
four relevant questions in this analysis: 

• do EPS trees produce the same range of a c values 
for a given mass as halo trees extracted from 
simulations? 

• do EPS trees produce the same trend of a c with 
M , which could then explain the trend of c v i r with 
M? 

• is a c (M) constant with redshift for EPS trees, 
which could then explain the c v i T (z) trend? 
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Fig. 9. — Probability distributions of a c and c v i r when other pa- 
rameters are fixed. Top left: the distribution of c v j r for a given 
value of a c ; the scatter is less than ~ 20% (with no correction for fit 
errors in e v j r or a c ), and less than ~ 15% when corrected for errors 
in c v ; r . Top right: the distribution of a c for halos in the mass range 
1.5 — 2.5 X 1O 12 /i _1 M . Bottom left: the distribution of c vir for 
halos in the same mass range. Bottom right: the distribution of 
c v i r for the same mass range, as predicted from a c via Eq. 7. The 
mean and scatter of each distribution is listed at the top of the panel 
and is virtually equivalent between the modeled c v ; r (bottom right) 
and the measured simulation values (bottom left). 



Fig. 10. — Comparison of the formation epoch a c as measured us- 
ing the trajectory to z = 1 and the trajectory after 2 = 1. Either of 
these values individually is well correlated with the standard value 
using the whole history (i.e., on average a c stays constant for a given 
trajectory), but they are not clearly correlated with each other, in- 
dicating that a halo's past does not predict its future. Note that 
the vertical scatter is larger than the horizontal scatter, reflecting 
the fact that the second half of the trajectory does not constrain the 
value of a c very well because this is the shallow part of the function. 



• does the scatter in a c for a given mass in EPS trees 
account for the scatter in the measured c v i r (Af) 
relation? 

8.1. a c vs M vir and the EPS Offset 

In order to answer these questions, we first generate 
a random ensemble of mass accretion histories. This is 
done by coupling the extended Press-Schechter formalism, 
which predicts the probability of accreting a given mass 
in a given time, with a method for generating Monte- 
Carlo merger trees. Here we use the scheme introduced 
by SK99, with an additional modification proposed by 
Bullock, Kravtsov, & Weinberg (2000). For completeness, 
we outline the fundamental aspects of the method in Ap- 
pendix C. 

We consider the list of all halos found in the sim- 
ulation with masses larger than the minimum-fit mass 
2 x 10 u /i -1 M Q (~ 3000 halos). We then generate ten 
Monte-Carlo realizations of mass accretion trajectories for 
each of these halos, based solely on its z — mass, keep- 
ing track of the growth of the most massive progenitor as 
a function of time. We start at z — and trace histories 
back to z — 7.2. Using this most massive progenitor tra- 
jectory, a is calculated by fitting the trajectory to Eq. 3. 
The value of a c is then defined by Eq. 4, with 5 = 2. 

Figure 11 shows the distribution of a c values found for 
different given mass ranges, for both the simulated halo 
trajectories and the EPS trajectories. We see that the 
distribution of a c values from EPS trees is systematically 



offset toward later formation times for all mass ranges, and 
also appear to be slightly broader. A similar discrepancy 
was found by SLKD, although they did not investigate 
the same quantity — they found the average mass of the 
largest progenitor to be larger in the EPS trees than in 
the simulated trees at low redshift, and smaller at high 
redshift. This implies later formation times for the sim- 
ulated halos, as seen here. This discrepancy can be un- 
derstood in terms of comparing the conditional progenitor 
mass functions; SLKD showed that EPS over-predicts it 
for low masses and under-predicts it for high masses, with 
the mass scale of the crossover decreasing with increas- 
ing redshift. This discrepancy in formation epochs seems 
to directly reflect the more well-known finding that PS 
over-predicts the number density of halos below M* and 
under-predicts the number density above M* (e.g., Gross 
et al. 1998; Sheth & Tormen 1999). Although only one 
cosmological model (ACDM) has been investigated here, 
the disagreement of the halo mass function with Press- 
Schechter is generic and we expect that the formation time 
discrepancy would be seen to some extent in any CDM 
cosmogony. 

Although we find a discrepancy in the median a c val- 
ues as a function of mass, we find that they are offset by 
a constant multiplicative factor. Figure 12 shows the me- 
dian and 68% scatter for a c (M), from the EPS trajectories 
and the simulated trajectories, for z — halos. The con- 
stant offset is in the sense that the characteristic formation 
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epochs derived from the EPS trajectories are roughly 25% 
larger than those measured from the simulation trees. The 
scatter is quite comparable, though it is slightly larger for 
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Fig. 11. — Comparison of a c values in the simulation with those 
derived from EPS trees, for various mass ranges in the different 
panels. The filled histogram represents simulated halos and the 
shaded histogram represents ten realizations of EPS trajectories for 
the same set of masses. Abnormal halos (whose trajectories end pre- 
maturely) in the simulations are counted in the left-most bin; the 
number of these cases is small, and negligible above 10 12 h~ 'Mg . 

the EPS-derived values. 

If we assume that the a c values found using EPS tra- 
jectories are correct except for this constant 25% offset, 
due to the known discrepancies in PS theory, then these 
values can be used in combination with Eq. 7 to estimate 
c vir for each halo, using its EPS-derived trajectory. Figure 
13 shows the measured c v i r (^Wvir) for simulated halos at 
z = 0, compared with the c V i r value obtained using this 
semi-analytic method. The constant c\ has been shifted 
by a factor of 1.25, but otherwise we are able to match 
both the mass trend and the scatter of the measured c V i r 
values in the simulation with this simple model. 

8.2. Redshift Dependence in EPS 

We can also test whether the redshift dependence of the 
EPS trajectories is consistent with our model, and with 
the c v j r trend seen in the simulation. B01 showed that for 
a given mass, c v - lr oc a, implying, for our model, that (on 
average) a c must be uniquely set by the mass, indepen- 
dent of redshift. Figure 14 shows the a c (M) relation at 
rcdshifts z = 0, 1,2, and 3, using one realization of EPS 
trajectories (again, with the mass weighting of the simu- 
lated halos). There is no discernible change in the median 
or scatter over the entire mass range examined. Using Eq. 
7 to translate these values into predicted c v - 1T values will 
thus result in exactly the trend found by B01, i.e., that 
c v ir oc a, with the correct mass trend and scatter for all 
redshifts. 

Given the fact that the a c values are robust to changes 
in z Q , it might be considered somewhat puzzling that 



the a c (M) relation remains constant with redshift, since 
clearly the halo masses are significantly lower at high red- 
shift. Naively, one might think that since the masses shift, 
the average halo of a given mass would have a higher a c 
value in the past. However, for a given mass halo at z = 0, 
those halos which have the latest formation times will have 
the smallest masses at high redshift, and those which form 
earliest will have the highest mass at high redshift. When 
halos are combined with the proper mass weighting, they 
thus conspire to keep the same relation regardless of red- 
shift. Figure 15 demonstrates this in detail. In the upper 
panel, the a c (M) trend is shown at z — for halos in two 
distinct mass ranges. In the lower panel, these same halos 
are shown at z = 3: the a c values for each halo remain 
essentially unchanged, but they are now plotted against 
their z = 3 masses. Halos in any given mass range at 
z = 3 consist of a combination number of low a c halos 
that will have low z = masses and high a c halos that 
will have high z = masses. These add in a manner that 
maintains the shape and normalization of a c (M). The fact 
that a c {M) is constant with redshift is an indication that 
the model proposed by B01 is a reasonable one — i.e., that 
(on average) the formation time of a halo is set only by 
its mass, and can be related to the time that mass was a 
fixed fraction of M* . 

In summary, the formation times of halos in EPS (us- 
ing an improved version of the SK99 method to gener- 
ate merger trees) are somewhat later than those found in 
the ART simulation. However, if we measure a c values 
for these semi-analytic trajectories (using Eq. 5), multi- 
ply these values by a constant factor (0.8), and translate 
these formation epochs to concentrations using the rela- 
tion Cvir = cia /a c , we are able, to a good approximation, 
to match the scatter, mass, and redshift dependence of 
the c v i r values found for simulated halos. This method 
can be used in semi-analytic models to estimate halo con- 
centrations that are both based on their individual mass 
growth histories and have the correct distribution at every 
redshift. 

The only disadvantage of the above method is that it 
requires generating full mass accretion histories for a sam- 
ple of halos. Since Eq. 5 is a one-parameter model, one 
might be tempted to directly calculate the more prevalent 
formation redshift Zfo.5, when the most massive progenitor 
mass was half of Mo, and whose probability distribution 
can be calculated directly from EPS without generating 
merger trees, and then translate this value into z c using 
Eq. 5 in order to derive a concentration. However, while 
this method predicts the mean values relatively well, in 
fact there is substantially more scatter in values of Zfo.5 
than in values of z c . This is due to the fact that indi- 
vidual trajectories are noisy, and at any one point in the 
trajectory, there is likely to be significantly more scatter 
than in the shape of the halo as a whole. 



9. SCATTER IN THE c v i r (M v ir) RELATION 

9.1. Evaluating a Corrected Scatter 

The scatter in the concentration parameter has been es- 
timated by B01 and Jing (2000), and it may have a num- 
ber of important observational implications. Jing (2000) 
found a scatter of 0.08 — 0.1 in log 10 c v ir, while B01 de- 
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Fig. 12. — Comparison of the median a c (M) relation with that 
derived from EPS trees. The triangles correspond to the trees from 
the ART simulation. The open circles represent ten Monte-Carlo 
realizations of EPS merger trees for each simulated halo mass. The 
EPS-derived a c values are offset about 25% higher than those mea- 
sured in the simulation. The filled circles are the EPS results shifted 
down by 20%. In each case the error bars are Poisson based on the 
number in the mass bin. The outer lines represent the 68% scatter 
about the median relation. The dark solid line in the middle repre- 
sents the model of B01 with F=0.015, and 35% scatter about that 
model. 
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Fig. 13. — Comparison of the median c v i r (M) relation for sim- 
ulated halos (solid lines) with those obtained using Eq. 7 and the 
a c values derived from EPS trajectories (dashed lines). As in Fig. 
12, the EPS etc values have been offset 20% lower. Ten Monte-Carlo 
merger tree realizations are generated for each simulated halo mass. 
Some selected realizations look virtually identical to the distribu- 
tion for simulated halos. The outer pair of thin lines represent the 
corresponding 68% scatter and error bars are Poisson based on the 
number of halos in each bin. 



rived a somewhat larger scatter of A log 1 



0.14. 



A revised scatter estimate from our improved halo cata- 
logs, which also relies on the mass trajectories, is presented 
here. We also discuss how the scatter in c v i r is affected by 
the merging history of halos. 

B01 devised a method for correcting the scatter esti- 
mate for errors in the fits, and actually plotted this cor- 
rected scatter in Figure 4. Our new halo catalogs have 
significantly fewer halos with large fit errors, and thus a 
smaller uncorrected scatter, but encouragingly, when we 
use the method of B01 for correcting this for fit errors and 
Poisson errors, we get almost identical results (shown in 
Figure 16). This method involves doing 500 Monte Carlo 
realizations of each halo, in which their c V i r value is cho- 
sen from a one-sided Gaussian deviate with a standard 
deviation error on the measured c V i r value. This value is 
then added or subtracted to the measured value depend- 
ing on whether the measured value was below or above the 
median value for that mass. Poisson errors due to finite 
statistics are then subtracted in quadrature to obtain an 
estimate for the intrinsic scatter — which we also find to 
be A(logc vir ) m 0.14. 

As mentioned previously, there may be halos which for 
a short period of time are not well-fit by an NFW density 
profile; in many cases this occurs when a halo is under- 
going merging or disruption. In most cases, this results 



in an artificially low value of c V i r , which usually persists 
only for one timestep. In order to remove the effects of 
badly estimated concentration parameters which persist 
only for a short time, we use the merging history of a halo 
to identify those cases where the concentration has jumped 
significantly since the previous output time. The c V h-(M) 
relation excluding these halos (which make up ~ 8% of 
the total) is plotted in Figure 16, compared to the whole 
sample; the scatter for this sample is reduced from ~ 38% 
to ~ 31%. Note that the exclusion of these halos does not 
change the median value significantly. We regard the scat- 
ter estimate from excluding these halos as a lower limit. 
With this correction, our analysis is closer to the scatter 
estimate presented by Jing (2000). It is possible that some 
of the remaining discrepancy could be due to an underes- 
timate of fit errors; the total scatter for the very mas- 
sive halos in our analysis is slightly smaller. However, it 
should also be pointed out that Jing (2000) only consid- 
ered relaxed halos, and our analysis contains halos with 
a full range of properties, including those that have been 
recently disrupted — and these add significantly to the 
scatter. 

9.2. Scatter and Halo Merger History 

As discussed in §5, halos with recent major mergers dis- 
play the same trend between concentration and a c (Eq. 
3), but with somewhat more scatter. We show in Fig- 

The 



5 The value of 0.18 reported in the abstract of B01 was actually the directly measured scatter, without any correction for fit errors 
corrected value of 0.14 is illustrated in B01 Figure 4. This footnote is aimed at correcting this oversight in that paper. 
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Fig. 14. — Robustness of a c (M) to the redshift of "observation" 
in one realization of an EPS tree. The 4 different sets of curves and 
symbols refer to 4 different redshifts. The thick lines represent the 
median a c values and the pairs of outer thin lines represent 68% scat- 
ter. Error bars represent Poisson error based on the number in each 
mass bin. As predicted by Eq. 3, the median a c value is uniquely 
determined by the halo mass, independently of the redshift. 
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Fig. 15. — An analysis of the robustness of a c (M) to redshift, 
in EPS trajectories. Top: a c (M) at z = 0, broken into two mass 
ranges. Bottom: a c (M) at z = 3 for the same mass ranges deter- 
mined at z = (dark thick: high mass; light thin: low mass). The 
ratio M(z=3)/M(z=0) is related to a c , such that for a given z = 
mass, halos with high a c will have lower 2 = 3 masses. The dashed 
line shows the sum of the mass bins, which follows a similar trend 
with mass as the 2 = halos. The apparent upturn in the left-most 
bin of the dashed curve is due to the exclusion of halos less massive 
than 4 X 10 /l — 1 M at z = 0, which contribute lower a c values. 



tire 17 that these halos also follow the same basic trend 
with mass. This figure compares this trend and scatter in 
c v i r (Af) for all halos with that seen in halos that have not 
had a major merger since z — 2. The scatter is reduced 
from about 31(38)% for all halos to - 26(28)% for those 
halos without recent major mergers, where the first listed 
estimate is when halos with large c vlT jumps are excluded, 
and the number in parenthesis includes all halos. (The 
implied correction for halos without recent mergers is not 
large since most of the halos with jumps arc in the process 
of merging.) The sample of halos which have had a major 
merger since z = 2 do not have reduced scatter compared 
to the whole sample. Since halos with recent major merg- 
ers have later formation times on average, when they are 
excluded from the sample, the remaining halos have higher 
c V ir values by a factor of about ~ 10% compared to the 
complete sample; this can also be seen in Figure 17. 

The amount of scatter in the concentration parameter 
for a given mass is of particular interest because of its pos- 
sible implications for scatter in the Tully-Fisher relation 
(see, e.g., B01; van den Bosch 2000), and there is particu- 
lar interest in the amount of scatter in the concentrations 
of spiral galaxies. In many scenarios for galaxy formation, 
major mergers destroy disks, and thus halos with recent 
major mergers are unlikely to host spiral galaxies. As de- 
scribed above, not counting these halos reduces the scat- 
ter of the whole sample and thus may reduce consequential 
scatter in Tully-Fisher to a level that can be matched with 
observations. Whether these remaining halos host spiral 



galaxies may be influenced by how much mass they have 
accreted since their last disruption; this could bring the 
concentrations of halos hosting spiral galaxies lower than 
the median shown here. It should also be noted that in this 
analysis we have only considered distinct halos. To get a 
full estimate of scatter and normalization for galactic halos 
we would have to include subhalos in the analysis, which 
B01 have shown have somewhat larger scatter in c v i r - We 
defer a full analysis of these issues to a later work. 

10. DISCUSSION AND CONCLUSIONS 

Making use of a large sample of dark halos simulated in a 
cosmological volume, we have studied the relation between 
mass accretion history and the density concentration of 
halos. Remarkably, halo mass growth curves (normalized 
to the final halo mass) can be accurately described by a 
one parameter function, in which mass accretion occurs 
rapidly at early times, and slows at late times. The char- 
acteristic "formation" time, a c , defined as the time when 
the log-mass infall rate drops below a fixed value, fully de- 
fines the trajectory. We find that the value of a c for each 
halo trajectory is independent of the epoch at which the 
halo is observed, a a . In addition, the average value of a c 
for halos of fixed mass is independent of redshift. 

The NFW halo concentration parameter, c V i r , which, 
in combination with the halo mass, uniquely sets the 
shape of the density profile, is tightly related to a c via 
c v ir = c\a /a c . The central density of a halo at fixed ra- 
dius grows rapidly when the mass accretion rate is high, 
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Fig. 16. — Scatter in the c vlr -M v i T relation for halos at z = 0. The 
thick lines represent all the halos. The thin lines represent all halos 
whose concentration has not jumped by more than a factor of two 
(in either direction) since the previous output time (at z = 0.01). In 
each case, the solid lines represent the median value of c v j r , plotted 
with Poisson error bars, and the dashed lines represent the scatter 
corrected for errors in the individual profile fits and for Poisson scat- 
ter in the bins. In the mass range where the results are most reliable 
(~ 1 X 10 12 /i _1 M e ), the value of the corrected scatter in these two 
samples is A(logc v j r ) 0.14 and 0.12, respectively. 
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Fig. 17. — Type dependence of the c v - lr -M v i T relation for halos at 
z = 0. The thick lines represent halos which have not had a major 
merger since 2 = 2; the thin lines represent all types. In each case 
halos with big 'jumps' in c v ; r are excluded, as described in the pre- 
vious caption. In each case, solid lines represent the median value 
of c v i r , plotted with Poisson error bars, and dashed lines represent 
the intrinsic corrected scatter. For halos without recent mergers, the 
scatter is reduced to A(logc v j r ) fs 0.10. 



and approaches a constant value as the mass accretion 
slows. This result is consistent with a picture in which 
the mass accretion rate determines how far accreted mass 
makes it into the center of a halo: for high mass accretion 
rates, accreted material makes it far into the center of the 
halo, but as it slows new material builds up on the outside. 
Thus central densities of halos asymptote to a value which 
is proportional to the density of the universe at the time 
when the mass accretion rate slows. In late forming halos 
this process is delayed and thus the final central density is 
lower. To demonstrate the model, in Figure 18 we plot the 
evolution of several variables for an early and late form- 
ing halo: the mass, concentration, log-slope of the mass 
accretion rate, scale radius R s , and the density within a 
fixed radius. Each parameter is calculated analytically; 
this is done using Eq. 1, which specifies the density pro- 
file, Eq. 5, which specifies the mass accretion history, and 
Eq. 7, which specifies the relationship between them. 

We showed that scatter in c V i r for a given mass can be 
explained almost exclusively by scatter in a c for halos of 
that mass. Thus this model, based on the c v i r -a c cor- 
relation, captures successfully the main properties of the 
concentration parameter, including its mass dependence, 
redshift trend, and the scatter about these relations. 

The formation times derived from semi-analytic real- 
izations of EPS mass accretion histories were compared 
with those measured in the N-body simulation, and found 
to be systematically larger by 25%. This offset reflects 
known inaccuracies in the Press-Schechter approximation. 



By adjusting the formation times for this offset, we have 
successfully used the c v - lr -a c correlation to reproduce the 
mass and redshift dependence of c V i r , and the scatter about 
these relations, using mass accretion histories derived from 
EPS. This technique is likely to be very useful for inclusion 
in semi-analytic models of galaxy formation, which so far 
have at best drawn c vlI randomly from an assumed global 
probability distribution, with no dependence on the halo's 
history. 

We have presented an estimate for the scatter in c V i r , 
with preliminary hints for how this scatter may depend on 
galaxy type. For halos of a given mass that have not had 
a major merger, we estimate that the scatter is as low as 
A(logc V i r ) w 0.1. This scatter is roughly consistent with 
the scatter observed in the Tully-Fisher relation. However, 
we have not included subhalos in this estimate, so our re- 
sults only apply to isolated field galaxies; including sub- 
halos may increase the scatter. It is also interesting that 
such halos (those without recent major mergers, which 
could conceivably host disk galaxies) are also somewhat 
more concentrated, which is perhaps contrary to naive ex- 
pectations. 

It is worth pointing out that this is not the only attempt 
to connect the full mass accretion history of halos to their 
density structure. For example, Avila-Reese, Firmani, & 
Hernandez (1998), building on the work of Zaroubi & Hoff- 
man (1993), used an analytic approach based on shell-by- 
shell mass accretion to connect halo structure with mass 
accretion (see also Ryden & Gunn 1987). Although their 
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Fig. 18. — Evolution of mass and structure of two halos. From 
top to bottom, we plot the evolution of the mass M V \ T , the concen- 
tration c v ; r , the log-rate of the mass accretion history dlogM/dloga, 
the scale radius R B , and the density of the halo within 10 kpc. For 
each parameter the evolution of two halos with the same final mass 
(5 X 10 12 ) but different values of a c (0.25 and 0.5) is shown. In the 
plot of the mass accretion rate, the value of S used to define a c is 
marked with a thin dotted line. 

results differ in detail from what we have found using N- 
body simulations, the general trends associated with early 
and late formation seem to agree. Our results may serve 
as a useful benchmark against which specific analytic and 
semi-analytic models of halo structure formation can be 



tested. 

The correlation we have found between individual halo 
assembly histories and their density profiles is likely to 
have important consequences for a large range of observ- 
able galaxy properties. For example, halo density profiles 
directly affect galaxy rotation curves, and are likely to play 
an important role in determining galaxy shapes, gas infall 
and star formation rates. We thus expect that the inclu- 
sion of such a correlation in the context of semi-analytic 
models which track galaxy formation and evolution with 
simple recipes will affect a number of predictions of these 
models, and thus will provide a significantly more realistic 
theoretical framework for understanding galaxy popula- 
tions, the origin of galaxy type and the variation in galaxy 
properties. 
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APPENDIX 

FINDING AND FITTING HALOS 

The details of our procedure for finding and fitting halos are as follows: 

1. We construct density field values by a Cloud-in-Cell (CIC) process (Hockney & Eastwood 1981) on the largest 
grid of the simulation AL, and rank the particles according to their local density as determined on this grid. We 
then search for the possible halo centers, using two sets of smoothing spheres; one, with a small radius, r sp i, in 
order to locate the centers of tight, small clumps; and the other, with a larger radius, r sp 2, in order to locate the 
centers of halos with diffuse cores. The larger radius, r sp 2, is set equal to i?™" 1 , the virial radius of a halo of mass 
M™™. The smaller radius is set to r sp i = /A™ 111 , a rough approximation to the radius within which our 

smallest halo would be expected to contain ~ 5 particles. 

For each set of spheres, we take from the ranked list the particle with the highest local density and place a sphere 
about its location. A second sphere is placed about the next particle in the list not contained in the first sphere. 
The process is continued until all of particles are contained within at least one sphere. Because we are only 
interested in centers of halos more massive than M™ r m , we discard each sphere with fewer than a set number of 
particles. The minimum number of particles required for a kept sphere is determined separately for each radius. 

For the r sp i spheres, we use the following conservative halo density profile: 

(where C is determined by fixing the minimum halo mass to be Af™ r m ), in order to estimate the minimum number 
of particles within r sp i: 

Nspl = l + 6[W% sp )V2-l]- (A2) 
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For the z = output of the 60/i _1 Mpc simulation we analyze, A sp i = 3 (rounding to the next lowest integer). 
Spheres of size r sp i with fewer than JV sp i particles are discarded. Similarly, all of the r sp 2 spheres containing fewer 
than A sp 2 = A™ ln particles are discarded. 

The final list of candidate halo centers is made up of all of the (small) r spl spheres, together with each of the r sp2 
spheres that does not contain an r sp i sphere. 

2. For each sphere of radius r sp = r sp i or r sp2 , whichever applies, we use the particle distribution within the sphere 
to find the center of mass and iterate until convergence. We repeat the procedure using a smaller radius, r = r^, 
where = r sp /2i . We continue this method until n = tl, where is defined by the criterion ri > 2/ rcs > r^+i, 
or until reduction leads to a sphere with fewer than jV sp i particles. 

3. We unify the spheres whose centers are within of each other. The unification is performed by making a density 
weighted guess for a common center of mass, and then iterating to find a center of mass for the unified object by 
counting particles. The size of sphere used to determine the center of mass is the smallest radius that will allow 
the new sphere to entirely contain both candidate halo spheres. 

4. For each candidate halo center we step out in radial shells of 3f res , counting enclosed particles, in order to find 
the outer radius of the halo: Rh = min(i? v j r , R t ). The radius i? v ; r is the virial radius, and Rt is a "truncation" 
radius, defined as the radius (< i? v ir) in which a rise in (spherical) density is detected ((ilogp/dlogr > 0). This is 
our method for estimating when a different halo starts to overlap with the current halo and is important for halos 
in crowded regions. We estimate the significance of a measured upturn using the Poisson noise associated with 
the number of particles in the radial bins considered. Only if the signal to noise of the upturn is larger than a^ t 
do we define a truncation radius. The value of o\R t is a free parameter. We use o\R t = 5. 6 

5. Among the halo candidates for which we have found an i? v ir> we discard those with A v ; r < N™ ln , where A v ; r is the 
number of particles within i? v ir- Among the halo candidates for which we have found a rise in spherical density, 
we discard those which contain less than N^ n particles, where N^ n = N™ ln if R t > i?™J. n , otherwise 

Ng* = N?*(Jh). (A3) 



vir 



The above constraint follows from an extrapolation of the minimum mass halo using an isothermal profile 
p(r) oc 1/r 2 . 

6. For halos with more than -/V™g" particles, we model the density profile of each halo using the NFW form (Eq 1) and 
determine the best fit i? s and p s values, which determine i? v ir and M v i r . The fitting procedure uses logarithmically 
spaced radial bins from max(2/ ros , 0.02 x niin(_R vir , i? t )) out to R h . If any bins are empty we decrease the number 
of bins by one until this is no longer the case. If the number of bins is reduced below three we discard the halo as 
a local perturbation. 

The fit takes into account the Poisson error in each bin due to the finite number of particles, and we obtain 
errors on the fit parameters (<r# s and a Ps ) using the covariance matrix in the fit routine. The errors on the fit 
parameters can be translated easily into errors for R v [ r {o~R wir ) and the estimated NFW mass of each halo, M v i r 
(cm)- In some cases, the fit does not converge. When this occurs, we mark the halo as a non-fit. This occurrence 
is rare for distinct halos, but is more common for subhalos (see below). This may reflect a tendency for subhalos 
defined with the current merging criteria to be poorly described by an NFW form, perhaps as a result of frequent 
interactions or close neighbors. 

7. We unify halos with centers that overlap by -Rcombine- For a given pair of halos with virial radii i? v ir,i, ^vir,2 5 
wc define this combination radius to be -Rcombine = min(_R V i r .i/2, i? v ir,i /2) . If either of the halos does not have 
a fit i?vir, we use the halo radius Rh in place of i? v j r . Our criterion is met if two (or more) halo centers are 
within i?combino of each other while at the same time having velocities which allow them to be bound to the 
common system. If such a case occurs, then along with the individual halo NFW fits, we fit another NFW profile 
about the common center of mass of the two combined halos and decide whether the candidate-unitcd-halos are 
bound/unbound to the common NFW fit using the radial escape velocity determined using the common NFW 
profile (see below). If both halos are bound we combine the two halos into one, and keep the common fit for the 
characteristic parameters. If at least one is not bound, we do not combine the halos. 

8. For each halo, we remove all unbound particles before we obtain the final fits. We loop over all particles within 
the halo and declare a particle at a distance r from the center of a halo to be unbound if its velocity relative to 
the center of mass velocity of the halo obeys v > \/2 |3>(r)| . For halos that have fits, we use the radial potential 

6 The choice was motivated by several tests using mock catalogs of halos in clusters designed to determine how varying <jR t affects our ability to 
fit the density profiles of subhalos. Although our results were not strongly dependent on this choice, wc did obtain the best fits using on t = 5. 



18 



WECHSLER ct al. 



for an NFW density profile 7 : 

log(l + x) 



$nfwM = -4nGp s R s 



(A4) 



2 _ 

X 

otherwise, we use the radial potential for a singular isothermal sphere with the same mass. 

After removal, we construct a new density profile (and find new NFW fit parameters if N p > A™g"). The 
procedure is repeated until the number of unbound particles becomes < 1% of the bound particles or until the 
total number of particles within the halo falls below N™ in . 

9. For each halo in the final catalog, we determine its NFW fit if N p > N™^, and record its fit parameters and their 
errors. We also measure and record its spin parameter, A, and the maximum of its circular velocity curve, V max . 

As described above, the halo catalog we have developed includes an arbitrary number of levels of substructure within 
halos. The full catalog with substructure should then include all halos directly around galaxies, above the relevant mass 
resolution, and thus is useful for a number of direct comparisons with observations. However, for many purposes, a halo 
catalog including only "distinct" halos, i.e., halos which are not subhalos of any larger halo, is sufficient and introduces 
significantly less complication. For this reason we have culled the full catalog into a smaller catalog that does not include 
subhalos within halos; this is the catalog analyzed here. 

Since there are multiple levels of substructure, the details depend slightly on the algorithm chosen. We take the 
maximum circular velocity to be the most reliable measure of the halo's size, since it is a measured quantity and doesn't 
rely on a fit, and can be defined equivalcntly for all halos. All halos whose centers lie within the virial radius of a larger 
halo are then designated as subhalos. Note that a halo that lies within the virial radius of subhalos is only removed if it 
itself is classified as a subhalo of a distinct halo. 



CORRECTING MASSES 

Our procedure of fitting density profiles is intended to give the best estimate possible of a halo's virial mass; however, 
it is subject to large errors when there are a small number of particles or especially when the halo is undergoing merging 
or disruption and is far from being a relaxed, spherical object. These uncertainties are taken into account in the fit errors, 
but for many purposes it is essential to have the best estimate possible of the halo's mass at each output time. Especially 
for consideration of the evolution of individual halo mass trajectories it would be useful to eliminate large jumps in the 
trajectories which are due to the above-mentioned irregularities and not to real changes in a halo's mass. In order to do 
this, for each halo we compare three masses: M v ; r , as measured from the NFW fit, Mh, the measured mass within Rh, 
and M tra j, which designates the mass interpolated between the most massive progenitor in the previous output time and 
the offspring halo in the subsequent output time (assuming they both exist). In most cases, M = Af v ; r (if M v ; r exists, 
otherwise M — Mh); it is only changed if this mass seems clearly inconsistent with the other mass estimates and does 
not seem reasonable. For most halos, the error on M v j r is small and M v ; r ~ Mh; in these cases M always equals M v ; r . 
However, if one of these is not the case, M v -„ is used if it is close to M tra j and otherwise Mh is used if it is close to M tra j. 
If neither seems consistent with the halo's trajectory, we use M = mcdian(M v ; r , Mh, M tra j). The details of the procedure 
are slightly more complicated, depending on the error on M v ; r , and we direct the reader to Wechsler (2001) for a complete 
description. 



GENERATING MERGER TREES 



For completeness, we outline here the fundamental aspects of EPS and our method for generating merger trees. Lacey 
& Cole (1993, hereafter LC93) introduced a method for calculating the probability that a halo of mass M accretes a given 
mass in a given time. Let S(M) = tr 2 (M) be the linear density variance on the mass scale M and w(t) = S c (t) be the 
linear density for collapsing structures at time t (see, e.g., White 1996). Given a halo of mass M at some time t , the 
probability that it accretes a mass AM in a time At is then 



(Aw) 2 
2AS 



dAS, (CI) 



where AS — S(M) — S(M + AM) and Aw = w(t) — w(t\ + At). In order to generate halo merging trees, one must 
implement this formula iteratively, with some algorithm for choosing progenitors. However, no method has been proposed 
that simultaneously matches the the conditional mass function and progenitor distribution of specified by Eq. CI exactly. 
We use the scheme suggested by Somerville & Kolatt (1999, SK99), which enforces mass conservation exactly but only 
reproduces the progenitor distribution of EPS approximately (for alternative techniques see Kauffmann et al. 1993 and 
LC93). In order to keep the trees finite, a minimum progenitor mass, M m is defined. Halo mass growth with AM < 
M m is treated as diffuse accretion. A key ingredient of this technique is that the timestep must be chosen such that 
Aw ,$ y/M m dS/dM in order to reproduce the expected conditional mass functions of extended Press-Schechter. Stepping 
back in time, the tree then provides a list of progenitors and their masses. In order to better match the analytic prediction 

7 Note that this potential is not necessarily the physical gravitational potential at the halo location. For a subhalo, for example, the host 
background potential is not included. 
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of the progenitor distribution, we apply an addition fix to the method of SK99, suggested by Bullock et al. (2000), which 
constrains the number of progenitors at any timestep to be close to the mean for that mass and redshift. 
n 
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